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Abstract 

Biogeography investigates spatial patterns of species distribution. Discontinuities in species distribution are identified as 
boundaries between biogeographic areas. Do these boundaries affect genetic connectivity? To address this question, a 
multifactorial hierarchical sampling design, across three of the major marine biogeographic boundaries in the central 
Mediterranean Sea (Ligurian-Tyrrhenian, Tyrrhenian-Ionian and Ionian-Adriatic) was carried out. Mitochondrial COI sequence 
polymorphism of seven species of Mediterranean benthic invertebrates was analysed. Two species showed significant 
genetic structure across the Tyrrhenian-Ionian boundary, as well as two other species across the Ionian Sea, a previously 
unknown phylogeographic barrier. The hypothesized barrier in the Ligurian-Tyrrhenian cannot be detected in the genetic 
structure of the investigated species. Connectivity patterns across species at distances up to 800 km apart confirmed that 
estimates of pelagic larval dispersal were poor predictors of the genetic structure. The detected genetic discontinuities 
seem more related to the effect of past historical events, though maintained by present day oceanographic processes. 
Multivariate statistical tools were used to test the consistency of the patterns across species, providing a conceptual 
framework for across-species barrier locations and strengths. Additional sequences retrieved from public databases 
supported our findings. Heterogeneity of phylogeographic patterns shown by the 7 investigated species is relevant to the 
understanding of the genetic diversity, and carry implications for conservation biology. 
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Introduction 

Through the study of species distribution as well as patterns of 
genetic structuring, a link between biogeography and phylogeog- 
raphy can be explored, which aids in the identification of 
discontinuities in the spatial distribution and genetic diversity of 
species [1]. The level of connectivity between geographic areas is 
mainly determined by the interplay between the dispersal capacity 
of a species (including pre and post-settlement processes [2]), the 
hydrological regime [3] and the occurrence of geological/ 
topographical boundaries, which may reduce the dispersal 
potential and hence the gene flow [4]. 

Major barriers to gene flow between oceanic basins have 
previously been identified (for example, Point Conception in 
California, [5]; the Eastern Pacific barrier, [6]; the Indo-Pacific 
Barrier, [7]; the South-eastern Australian barrier, [8]). These 
barriers support the occurrence of genetic discontinuities between 
the world's biogeographical regions [9]. Significant genetic 
structuring has also been detected at smaller spatial scales in 
correspondence with the boundaries between biogeographic 
provinces [10] or sectors [11,12]. 

In the Central Mediterranean Sea nine biogeographic sectors 
have been proposed based on species' distributions [13], but their 
consistency with phylogeographic sectors has only been investi- 
gated in a few cases [14—16]. These studies investigated genetic 
differentiation between the Mediterranean basins (eastern, west- 



ern, Adriatic Sea) in a single species. However, data from a single 
species, with specific life history traits, ecology, ethology, etc, 
cannot be used to generalize phylogeographic barriers or 
limitations to gene flow across other species. In fact, studies on 
different target species conducted in the same geographical areas 
often reveal contrasting genetic structures and/or phylogeogra- 
phical patterns [17-19]. Comparative phylogeography has 
recently been suggested as a powerful tool to discern general 
patterns over species-specific traits [20,21]. It has proven effective 
in identifying across taxa barriers to gene flow in the Almeria- 
Oran Front [22], and along the north-eastern Pacific coast [23]. 
Currently no multispecies studies identifying barriers to gene flow 
or genetic discontinuities in the Mediterranean Sea exist. Given 
the practical implications of defining genetic units for the 
conservation and management of natural resources, a study of 
this nature is sorely needed [24] . 

In the present study, the occurrence of phylogeographic 
boundaries in the Central Mediterranean Sea, and their consis- 
tency across species was investigated. By means of a multifactorial 
hierarchical sampling design, three putative long-term barriers to 
gene flow were tested. These putative barriers were established 
based on known biogeographic boundaries [13]. Two replicated 
areas were sampled across each barrier, and within each area the 
genetic structure of seven common shallow water invertebrate 
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species was analysed using a portion of the COI mitochondrial 
gene. 

In addition to genetic differentiation measures, two types of 
multivariate analysis were used to test the consistency of the 
phylogeographical patterns across species. The results of this study 
provide a first contribution towards a comparative phylogeogra- 
phy in the Central Mediterranean Sea, with the aim of fully 
understanding the main past and present drivers of genetic 
diversity and differentiation. 

Materials and Methods 

Study areas and species 

Three major potential biogeographic/phylogeographic bound- 
aries, located in the Central Mediterranean Sea (Figure 1), and 
hereafter referred to as putative barriers (PB), were investigated. 
The first putative barrier (PB1, in Location 1) separates the 
Ligurian Sea and the North Tyrrhenian Sea, corresponding to a 
virtual line between the coast of Tuscany, and the islands of Elba 
and Corsica. The two sides of PB1 are characterized by divergent 
currents [25]. The second putative barrier (PB2, in Location 2) 
separates the Tyrrhenian and Ionian Seas, and corresponds to a 
line across the Strait of Messina and Sicily. Very strong tidal 
currents characterize the narrow Strait of Messina, flowing mainly 
from north to south [25]. The third putative barrier (PB3, in 
Location 3) separates the Ionian and Adriatic Seas, around the 
Apulian Cape of Santa Maria di Leuca. Strong currents with 
different patterns in opposite directions characterize the water 
circulation along this putative barrier [26]. 

At each side of the putative barriers two replicate sites 1 50 km 
apart from each other were defined (sites A, B, C, D, E, F). In each 




Figure 1. Map of the study area, locations defined, sampling 
areas and putative barriers analyzed. Location 1: A1. Castiglion- 
cello, A2. Calafuria, B1. Talamone, B2. Porto Ercole; Location 2: CI. 
Nicotera, C2. Vibo Marina, D1. Soverato, D2. Copanello; and Location 3: 
El. Gallipoli, E2. Porto Cesareo, F1. Porto Badisco, and F2. Otranto. PB1: 
putative barrier between Ligurian and Tyrrhenian Seas, PB2: putative 
barrier between Tyrrhenian and Ionian Seas and, PB3: putative barrier 
between Ionian and Adriatic Seas. Grey arrows represent the direction 
of the dominant currents in those areas. The bottom left image is taken 
from [65] and shows the most commonly acknowledged biogeographic 
sectors in the central Mediterranean region. 
doi:1 0.1 371 /journal.pone.01 01 1 35.g001 



site 2 sampling areas about 25 km apart were randomly selected, 
and named counter-clockwise starting from Location 1 (areas Al, 
A2; Bl, B2; CI, C2; Dl, D2; El, E2; Fl, F2; Figure 1). 

Seven invertebrate species were selected based on their 
prevalence in the central Mediterranean as well as to cover a 
variety of life history traits. The selected species belong to different 
taxa: Gastropoda (3), Polyplacophora (1), Ascidiacea (1); Crustacea 
(1), and Porifera (1). Data retrieved on life history traits and larval 
ecology are summarized in Figure 2. The longest Pelagic Larval 
Duration (PLD) estimates correspond to Balanus perforatus, Patella 
caerulea, and Halocynthia papulosa. Hexapkx trunculus is a brooding 
species and has no theoretical PLD. PLD estimates for Chondrosia 
rentformis, Chiton olivaceus, and Osilinus turbinatus are of about one 
week. 

In each area 30 individuals per species were collected, fixed in 
80% ethanol and maintained at 4°C prior to processing. Only the 
sponge was sampled as pieces rather than as whole organism. To 
minimize the possibility of clonality, samples were collected at least 
5 m apart [27]. No specific permits were required for the areas 
and the species studied. 

DNA extraction and COI amplification 

Total DNA was extracted from a piece of circa 5 g of tissue 
using a REDExtract-N-Amp kit (Sigma-Aldrich). The COI gene 
was chosen for analysis due to its intermediate rate of variation 
that provides information on genetic structure across a wide range 
of spatial and temporal scales [28] . Universal primers described in 
[29] were initially used for the amplification of a fragment of the 
COI gene. After the first sequences were retrieved, six pairs of 
species-specific primers were designed with the software PRIMER 
v. 3.0 (http://www.fokker.wi.mit.edu/primer3/input.htm) in or- 
der to increase amplification success and the quality of sequences 
obtained (Table SI). 

PCR amplification reactions were carried out in a 25 ul total 
volume with 5 u.1 of 5x GoTaq Flexi Buffer (Promega), 4 ul of 
MgCl 2 25 mM, 0.5 ul of dNTP mix 10 mM, 0.5 \x\ of each 
primer, 0.15 |J,1 of GoTaq DNA polymerase (Promega) and 2.5 |il 
of template DNA (from a 1:50 dilution of the extracted DNA). A 
hold at 94°C for 3 min was followed by 30 cycles of denaturation 
at 94°C for 45 s, annealing at a primer specific temperature for 
45 s and extension at 72°C for 90 s, finishing with a final extension 
at 72°C for 5 min on a thermal cycler (Applied Biosystems, 
GeneAmp 2700) (Table SI). PCR products were purified and 
sequenced by MACROGEN INC. 

All sequences were checked and edited manually on Chromas 
Pro (Technelysium, Pty, Ltd, Australia). Sequence alignments 
were performed using Mega v. 5 [30]. 

Data analysis 

Preliminary analysis showed no differentiation in any of the 
investigated species between the two replicated areas of each site. 
As a result, individuals from the two areas were pooled and 
hereafter data have been analysed considering only the sites (A, B, 
C, D, E, F; Table 1). We detected no genetic differentiation 
between Locations 1 and 3 of Balanus perforatus, 1500 km apart, so 
decided not to analyse Location 2 samples. 

a. Genetic diversity, differentiation and gene flow 
estimates 

Number of haplotypes (h), haplotype and nucleotide diversity 
(Hd and 7t, respectively), were calculated for each species by site, 
location, and for the whole study area, using DnaSP v. 4.50.3 [31]. 
Diversity values and relationships between the haplotypes were 
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visualized in unrooted haplotype networks, calculated by Median 
Joining with the software Network v. 4.6.1.1 [32]. Loops were 
solved when possible following criteria derived from coalescent 
theory [33]. 

The fixation index (Fst) between sites across each PB was used 
to quantify the differentiation due to genetic structure. Large scale 
genetic structuring across the whole study area was also assessed, 
as the genetic differentiation between all pairs of sites. Pairwise Fst 
values were calculated in Arlequin v 3.1 [34] with 10000 
permutations, and their significance corrected for multiple testing 
following [35]. 

Analysis of Molecular Variance (AMOVA) was performed to 
investigate sources of variability within and between the three 
locations studied. Significance was tested with 16000 permutations 
in Arlequin v 3.1 [34]. 

Isolation by distance across the whole study area was tested with 
a Mantel test as implemented in Genepop, following the Isolde 
procedure which compares a semi matrix of Fst/ ( 1 -Fst) with a semi 
matrix of the natural logarithm of the shortest oceanic distance in 
km. 

Migration rates between areas were estimated with the software 
MIGRATE-n v. 3.3.1 [36], using Markov chain parameters. 

COI sequences of the target species available in Genbank 
(Table S2) were aligned and represented in haplotype networks 
with the haplotypes obtained in this study, to provide some insight 
on large-scale genetic diversity of the species and patterns across 
other acknowledged barriers to gene flow. 

b. Estimation of effective population size variations 

The history of effective population size (Ne) was assessed by 
testing the fit of the mismatch distributions to the theoretical 
distribution in an expansion scenario [37] using Arlequin. Fu's and 
Tajima's neutrality tests were also carried out, as they are more 
sensitive to population expansions than the mismatch distributions 
[38]. These tests were performed for each site, for the whole data 
set, and for undifferentiated groups of sites identified by the 
previous population differentiation analysis. 

For those sites that were not significantly different, and did not 
deviate from the expansion model, the coalescence time X (i.e. time 
since the start of a population expansion) was estimated following 
the formula T = 2ut [37], where T is the date of growth or decline 
in units of mutational time and u = 2 uk, where k is the number of 
generations per year and [i is the mutation rate per nucleotide. 

For the same undifferentiated sites, past changes in effective 
population sizes and estimates of expansion time were also inferred 
by generating Bayesian Skyline Plots (BSPs) using the software 
BEAST v. 1.7.5 [39]. The best-fit models of nucleotide substitu- 
tion for each species data set were selected with jModeltest v. 0.1.1 
[40] and fed into BEAST. Mutation rate estimates used were: 
2.4% MY for molluscs [41], 3.1% MY for barnacles [42], and 
2.86% MY for ascidians [43]. Given the absence of molecular 
clock estimates for sponges, the BSP for this species was not 
calculated. Confidence intervals for Ne through time were 
calculated from the posterior probability distributions using 
TRACER [44]. 

c. Multispecies analysis 

To test for consistency of the phylogeographic patterns across 
species, two types of multivariate analysis were used. 

The multivariate Analysis of Similarity (ANOSIM), as imple- 
mented in Primer v.6 [45], was used to test the differences across 
each of the three putative barriers. With this analysis we aimed to 
provide an across-species analysis of haplotype frequencies across 
the putative barriers, reassuming most of the information on the 
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Figure 2. Pelagic larval dispersal estimates, spawning period, and other reproductive data for the species studied. Data retrieved 
from [58,66-74]. Roman numerals correspond to the months and black bars indicate the extent of the spawning season. * The accepted name of 
Osilinus turbinates is Phorcus turbinates. ** The accepted name of Balanus perforates is Perforates perforates (World Register of Marine Species, www. 
marinespecies.org). 
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genetic structure of all species, and minimizing the effect of 
between species sample size and/or genetic diversity differences. 
The frequencies of each haplotype of all seven species were the 
variables, and the Areas were the samples. In order to accomplish 
with the independence of the variables, the haplotype frequencies 
were calculated for each haplotype across all sites (instead of for 
each site across all haplotypes). We used Site and Location as crossed 
factors in order to test for the effect of each of the PB. A Bray 
Curtis similarity matrix between all sites was represented 
graphically using a Multidimensional Scaling plot (MDS). 
Apparent groupings across the whole study area were also 
statistically tested with ANOSFM. 

Discriminant analysis of principal components (DAPC), avail- 
able in the Adegenet package [46] for R (R core team 2012), was 
also performed. This technique is designed for multivariate genetic 
data (multi-locus genotypes). It maximizes the variation between 
groups by first performing a principal component analysis (PCA) 
on pre-defined groups or populations and then uses the PCA 
factors as variables for a discriminant analysis (DA), thus ensuring 
their independence. To convert our sequences into a multivariate 
dataset, each polymorphic position of each sequence was codified 
with a four-digit code, corresponding to the four nucleotides. Only 
the polymorphic sites were considered, creating a "pseudo- 
genotype" of each specimen. Polymorphic sites across species 
were counted as null alleles in order to test all species 
simultaneously. In this case, the number of variables corresponds 
to the total number of specimens and the samples correspond to 
the six sites. This analysis allowed us to take into account not only 
the haplotype frequencies differences, but also how different these 
haplotypes are, minimizing the effect of the species (although 
unavoidable due to the different levels of polymorphism, and 
hence the number of "pseudo-loci" included per species). 



Results 

Genetic diversity and connectivity patterns 

A total of 199 haplotypes were deposited in Genbank with 
accession numbers KF297371 - KF297570. The number of 
haplotypes discovered ranged between 5 in Chondrosia reniformis and 
48 in Osilinus turbinates. Patella caerulea and Balanus perforatus 
presented 36 and 43 haplotypes, respectively. Hexaplex trunculus, 
Chiton olivaceus, and Halocynthia papulosa showed lower number of 
haplotypes (Table 1). The highest values of nucleotide diversity 
were found in H. trunculus, and H. papulosa; the highest haplotype 
diversity values were found in B. perforatus and H. papulosa. C. 
reniformis and C. olivaceus showed the lowest diversity values 
(Table 1). H. trunculus and P. caerulea both showed the highest 
diversity values at Location 2. H. trunculus at Location 3 also 
showed high diversity values, whereas P. caerulea showed higher 
values at Location 1. 

Patella caerulea and Hexaplex trunculus showed two haplogroups in 
their networks. P. caerulea presented two common haplotypes 
separated by a single mutation and surrounded by several low 
frequency haplotypes, 26 of which were private (Figure 3a). In H. 
trunculus, nine mutation steps separated the haplogroups, which 
could therefore be considered as two divergent lineages. The most 
common lineage occurred across all sites, whereas the second 
lineage was found in sites D, E, and F (except for H_27, Figure 3b). 
Two other highly divergent haplotypes were found in site A 
(H_16) and site E (H_23). 

Three species presented star-like shaped haplotype networks 
with different degrees of ramification. The Balanus perforatus 
network (Figure 3g) showed a high number of low frequency 
haplotypes. Osilinus turhinatus also showed a profusely branched 
network, with several middle frequency haplotypes (H_l, H_ll, 
H_l 3, H_3 1 ; Figure 3c). The Chiton olivaceous network had a simple 
star shape: one very common haplotype surrounded by several 
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Table 2. Pairwise Fst values between sites for all seven species. 
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0.22265 


0 








D 0.67033 


0.66462 


0.53127 


0 






E 0.59243 


0.58192 


0.48037 


0.02646 


0 




F 0.37615 


0.34750 


0.26302 


0.14552 


0.03280 


0 


c. Osilinus turbinatus 


A 0 


B -0.00414 


0 










C 0.08221 


0.09829 


0 








D 0.01355 


0.04912 


0.04474 


0 






E 0.09587 


0.10261 


0.16009 


0.15478 


0 




F 0.09282 


0.10546 


0.10968 


0.08825 


0.03916 


0 


d. Chondrosia reniformis 


A 0 


B 0.0318 


0 










C -0.0293 


0.0274 


0 








D -0.0293 


0.0274 


0.0000 


0 






E 0.4632 


0.3062 


0.3838 


0.3838 


0 




F 0.4739 


0.3631 


0.4078 


0.4078 


-0.0509 


0 


e. Chiton olivaceus 


A 0 


B -0.0073 


0 












C -0.0136 


0.0000 


0 








D 0.0071 


0.0143 


-0.0167 


0 






E 0.0653 


0.0577 


0.0512 


-0.0172 


0 




F 0.1428 


0.1309 


0.1405 


0.0005 


0.0118 


0 


f. Halocynthia papulosa 


A 0 


B 0.0152 


0 










C 0.1212 


0.1414 


0 








D 0.2030 


0.2225 


0.0890 


0 






E 0.0903 


0.1705 


0.2325 


0,1024 


0 




F 0.0897 


0.1716 


0.2958 


0.2018 


-0.0336 


0 


g. Balanus perforatus 


A 0 


B -0.0127 


0 










C 


D - 


E -0.0206 


-0.0071 






0 




F -0.0188 


-0.0137 






-0.0201 


0 



Values in bold indicate significance after correction for multiple testing. 
doi:10.1371/journal.pone.0101135.t002 
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private and three low-frequency haplotypes, with a maximum 
distance of three mutation steps (Figure 3e). 

The five haplotypes of Chondrosia reniformis (one common 
haplotype, three privates, and a middle frequency haplotype in 
sites E and F, Figure 3d) were separated by a maximum of 2 
mutation steps. Halocynthia papulosa presented several closely 
related common haplotypes, and few private haplotypes (11 out 
of 21, Figure 3f). 

Fst values failed to detect significant genetic differentiation 
across the putative barriers PB1 and PB3 for any of the seven 
species studied (Table 2). Significant genetic differentiation (p< 
0.015 after FDR correction) was only detected across the PB2 for 
Patella caeruka and Hexaplex trunculus (Table 2a— b). Pairwise Fst 
values across the whole study area showed significant differences in 
P. caeruka and H. trunculus between northwestern (A, B and C) and 
southeastern sites (D, E and F). P. caeruka migration rate estimates 
were higher from the southeastern towards the northwestern sites 
(M N . >S = 96.8; M S _ >N = 795.0), whereas H. trunculus migration 
rate estimates were higher from northwestern towards southeast- 
ern sites (M N _ >S = 658.5; M S .> N = 25.4). The Mantel test on the 
geographic and genetic distances matrix was not significant 
(p = 0. 1 35 for P. caeruka and p = 0.082 for H. trunculus). 

In Osilinus turbinatus and Chondrosia reniformis significant pairwise 
Fst values were detected between the E and F sites (Location 3) 
and most other sites (Figure 3c-d), except sites D and F in O. 
turbinatus, and sites B and E in C. reniformis (Table 2c-d). Migration 
rate estimates were bidirectional but asymmetric, higher from sites 
A, B, C, D, towards sites E and F (0. turbinatus M 1+2 .> 3 = 760.0, 



A x 

S3 


o 


2D Stress: 0 


O 








Figure 4. Graphic representation of the multivariate analyses. 

a. MDS of the Bray-Curtis distance matrix of sites according to the 
distribution of haplotype frequencies across all species, b. Graphic 
result of DAPC. Colour codes are the same as in Figure 3. 
doi:1 0.1 371 /journal.pone.01 01 1 35.g004 



M 3 _>i + 2 = 85.2; C. reniformis M 1+2 . >3 = 656.1, M 3 _> 1+2 = 363.0). 
Mantel tests were not significant in either species (p = 0.136 and 
0.101, respectively). 

Chiton olivaceus showed significant genetic differentiation between 
sites A, B, C and site F. Site E was also significantly different from 
site A (Table 2e). A significant Mantel test (p = 0.014) suggests 
that this pattern of genetic differentiation could be attributed 
to isolation by distance. Migration estimates were higher 
between nearby than between distant locations (M 1 _ >2 = 142.0; 
MPB 1 .> 3 = 154.7; M 2 .>i = 353.6; M 2 _ >3 = 81 1.6; M 3 _>i = 315.4; 
M 3 _ >2 = 302.3). 

Halocynthia papulosa showed significant genetic differentiation 
between several pairs of sites (Table 2f,) with no relation to 
distance (Mantel test p = 0.158). AMOVA analysis on the three 
locations wasn't significant (p = 0.064; Table S3), however it 
attributed a high percentage of variation to the differences 
between them (Va% = 13.56). Migration rate estimates showed 
wide confidence intervals (data not shown). The acorn barnacle 
Balanus perforatum showed no genetic differentiation between sites 
across putative barriers (PB1 and PB3), or even across sites located 
more than 1000 km apart (Table 2g). 

Among the 32 sequences of Patella caeruka retrieved from 
Genbank (Table S2), 12 corresponded to the H_4 haplotype (7 
from the Western Mediterranean, 2 from Lampedusa, and 3 from 
west of the Almeria-Oran Front) and six sequences corresponded 
to the H_l haplotype (from Greece, Lampedusa, Tuscany, 
Valencia, and 2 from west AOF). The haplotypes H_ll, H_26 
and H_17 (Figure 3a, northwestern haplogroup) were found in 6 
individuals from the western Mediterranean region, and 4 
haplotypes were not found in our dataset. The 1 1 COI sequences 
of Hexaplex trunculus retrieved belonged to individuals from the 
South of Portugal (Table S2). Nine of them corresponded to the 
H_4 haplotype and the other two haplotypes were one nucleotide 
different from H_4. 

Among the 10 Osilinus turbinatus COI sequences retrieved (Table 
S2), 3 (2 from Sardinia and one from SE Spain) presented the H_4 
haplotype (Fig. 3c); 3 corresponded to H_9 (Turkey, Cyprus and 
Croatia); and 2 belonged to H_24 (SE Spain). Two individuals 
from Croatia and one from Cyprus presented private haplotypes. 

The single sequence from Chondrosia reniformis published in 
Genbank corresponded to the H_2 haplotype and belonged to an 
individual from Israel (Fig. 3d). The Chiton olivaceus sequence came 
from an individual from north-eastern Spain and presented the 
H_2 haplotype (Fig 3e). The six sequences of Halocynthia papulosa 
(Table S2) belonged to individuals from the north-western region 
of the Mediterranean. Two of these presented the H_12 
haplotype, while the other two presented the H_13 haplotype 
and the last two showed the H_18 haplotype (Fig. 3f). 

Estimation of population size variations 

In four species, mismatch distributions and neutrality tests 
indicated expansion in most sites, as well as in the whole study area 
(Figure SI; Table S4). The brooding gastropod Hexaplex trunculus 
did not show signs of expansion in most sites (Table S4b), probably 
due to the admixture of two genetic pools, visible in the mismatch 
distribution (Figure Sib) and previously detected as the divergent 
lineages in the haplotype network. The solitary ascidian Halocynthia 
papulosa and the sponge Chondrosia reniformis did not show signs of 
expansion (Table S4). 

Bayesian skyline plots (BSP) showed large confidence intervals 
in most cases (Figure S2), yet they roughly reflected the estimates 
of time since expansion calculated with Rogers & Harpending's 
formula. The northern haplogroup of Patella caeruka and the 
southern haplogroup of Hexaplex trunculus demonstrated stable 
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Table 3. Bray-Curtis distance matrix between sites according to haplotype frequency distribution across all species. 





A 


B 


54.93 










C 


41.95 


42.55 








D 


30.96 


32.69 


35.41 






E 


34.97 


32.12 


32.85 


41.38 




F 


31.76 


29.93 


30.51 


39.34 


64.12 



doi:1 0.1 371 /journal.pone.01 01 1 35.W03 



population sizes (Ne) through time, coinciding with estimates of 
times since expansion (43,600 and 68,200 years ago, respectively) 
largely preceding the Last Glacial Maximum (LGM: 24,000- 
27,000 years ago). The more recent lineages showed a slightly 
growing population size in the BSP (Figure S2), but still dating 
before the LGM: 29,700 years ago in the southern haplogroup of 
P. caerulea and 31,700 ya in the northern lineage of H. trmculus. 
The BSP of Balanus perforates, Orilinus turbinates, and Chiton olivaceus 
showed increasing Ne through time consistent with the strong 
signs of expansion detected. Estimates of time since expansion 
were in accordance with an old expansion for B. perforates 
(62,600 ya), and coincided with the LGM for O. turbinates 
(25,300 ya). Expansion estimates for C. olivaceous dated around 
32,600 ya, before the LGM, in contrast with the BSP, which 
showed increasing population size after the LGM. 

Multispecies analysis 

A matrix of 159 different variables corresponding to the 
frequencies of the haplotypes of all species (excluding Balanus 
perforates) across sampling areas was generated. The analysis of 
similarity (ANOSIM) indicated no significant effect of the factor 
site crossed by location, i.e. no multispecies significant differences 
across any of the three putative barriers. The MDS of the Bray- 
Curtis matrix of distances between all six sites (Table 3) showed 
higher similarity of sites A, B (Location 1) and C (Location 2). Site 
D (Location 2) was isolated, and sites E and F (Location 3) 
clustered together (Figure 4a). ANOSIM on these three groups 
showed significant differences between them (R= 1; p = 0.017). 

After transforming haplotypes into pseudo-genotypes, we 
obtained a matrix of 802 individuals with genotypes of 152 loci 
each (most of them zeros), corresponding to all polymorphic sites 
across all species (33 of Patella caerulea, 50 of Hexaplex truncutes, 8 of 
Halocynthia papulosa, 37 of Osilinus turbinates, 19 of Chiton olivaceus, 
and 4 of Chondrosia reniformis). The resulting DAPC showed a 
distribution of the samples (sites) similar to the MDS plot: sites A, 
B and C grouped together. Site D was separated from all the 
others, and sites E and F clustered together (Figure 4b). 

Discussion 

In the study area two zones of genetic discontinuity have been 
identified: across the biogeographic boundary between the 
Tyrrhenian and Ionian Seas, and between the western and 
eastern coasts of the Ionian Sea. 

The strongest genetic structure detected was across the 
Tyrrhenian-Ionian biogeographic boundary. The presence of this 
long-term boundary was reflected in the genetic structure of two 
species out of the seven studied. It probably corresponds to the 
boundary between biogeographic sectors 3 and 6 as defined by 
[13] (see Figure 1). The genetic structure across these sectors has 



been described by [47—49], but the location of the boundary 
remained vague. In this area, due to the complex topography of 
straits and islands, Pleistocene sea level fluctuations have left a 
signature in the genetic structure of many species [48], subjected 
to periodic divergences and reconnections of the populations on 
the two sides of the Messina Strait. Past history events might 
explain the origin of haplogroups or divergent lineages, whereas 
present day events, such as the current regimes, aid in maintaining 
and enhancing the differences due to limited or asymmetric 
connectivity between sites. 

Patella caerulea showed significant genetic structure across this 
Tyrrhenian-Ionian PB, reflected in the two haplogroups with 
different frequencies across each side of the boundary. In spite of 
the significant differentiation, a bidirectional but asymmetric 
migration across the barrier was estimated. Other limpet species 
showed comparable genetic structure in this area [50] . In addition, 
the sequences retrieved in Genbank gave further support to this 
structure: western and eastern Mediterranean samples clustered 
more often with haplotypes of the "northern haplogroup", and 
"southern haplogroup", respectively. It is noteworthy that 
individuals from the Atlantic side of the Almeria-Oran Front 
(AOF), the best-known barrier to gene flow in the Mediterranean 
[20], presented common haplotypes, suggesting this barrier had no 
large effect on P. caerulea genetic structure [51]. Expansion 
estimates indicate long-term persistence of this species and, in 
accordance with the older expansion estimates, an older origin of 
the northwestern haplogroup. This is in agreement with the 
widespread distribution and abundance of this species and with the 
more probable Atlantic origin of this genus [52]. 

The subtidal gastropod Hexaplex trunculus, a brooding species, 
showed the sharpest genetic structure across the Tyrrhenian- 
Ionian biogeographic boundary. Despite this fact, several haplo- 
types were common to sites located up to 1500 km apart, on both 
sides of the putative barrier, probably indicating high migration 
rates from the Ligurian-Tyrrhenian towards the Ionian-Adriatic 
populations, but almost none in the opposite direction. This sharp 
asymmetry has been described in other regions and species, and is 
suggested to reflect persistent current regimes affecting larval 
dispersal during the reproductive period [53]. Estimates of 
expansion of both lineages are consistent with the persistence of 
this species during the LGM, in more than one glacial refuge, 
followed by a unidirectional secondary contact [54]. Expansion 
estimates are older in the southern lineage, which, together with 
the higher diversity values in this area, indicates a major 
contribution of the Eastern Mediterranean genetic pool to the 
total genetic diversity of this species [41]. All sequences retrieved 
from Genbank belonged to individuals from South Portugal, so no 
further details on the structure of other Mediterranean sites could 
be inferred. Most individuals showed haplotypes of the northern 
haplogroup, reinforcing their Mediterranean origin, but also 



PLOS ONE | www.plosone.org 



9 



July 2014 | Volume 9 | Issue 7 | e101 135 



Genetic Structuring across Biogeographic Barriers 



suggesting no effect of the AOF on the long-term genetic structure 
of this species, as was the case for Patella caerulea. 

An unexpected area of genetic discontinuity was disclosed 
between the eastern and western Ionian Sea. Two species, 
presumably not affected by the Tyrrhenian-Ionian boundary 
previously described, showed significant genetic structure across 
the Ionian Sea, i.e. between the Adriatic-Ionian sites and the rest 
of the study area. The Ionian Sea has never been considered to 
host phylogeographic barriers, either related to historical events 
nor according to present day oceanographic processes. No records 
referring to genetic structure across this region were found in the 
literature, despite the fact that the genetic structures of some 
species support this finding [48], although no explicit tests to 
identify the barrier had been carried out. The location of the 
biogeographic boundary separating the eastern and western 
Mediterranean is controversial [55]. This newly detected phylo- 
geographic boundary provides new insight to understanding the 
biogeography of the area, although further data on the genetic 
structure of other species in the area would be needed. 

Osilinus turbinatus presented a high diversity of haplotypes and 
significant structure across the Ionian Sea, with bidirectional and 
symmetric migration estimates. Published sequences supported the 
detected structure and individuals from the eastern Mediterranean 
and the Adriatic Sea presented haplotypes related to those found 
exclusively in the Adriatic-Ionian sites. The estimates of expansion 
of this species coincide with the Last Glacial Maximum (LGM 
between 24,000 and 27,000 ya, [20]). 

The sponge Chondrosia reniformis, in spite of evincing low diversity 
values, also showed significant genetic structure across the Ionian 
Sea. A similar differentiation has been described in other 
widespread sponge species in the same geographical area [56]. 
The only sequence retrieved from Genbank belongs to an eastern 
Mediterranean individual and shares the more common haplo- 
type. Although further data from the Ionian Sea and the eastern 
Mediterranean would be needed to confirm our finding, this study 
suggests that for Osilinus turbinatus and C. reniformis the phylogeo- 
graphic barrier between eastern and western Mediterranean is 
located in the Mid Ionian Sea. 

Significant genetic structure was detected across the northern 
and southern Tyrrhenian Sea in Hexaplex trunculus and Osilinus 
turbinatus. In both cases significant Fst values were likely due to the 
occurrence of several private haplotypes, south Tyrrhenian 
haplotypes in the case of H. trunculus, and Ligurian-North 
Tyrrhenian private haplotypes in the case of 0. turbinatus. 

Halocynthia papulosa showed several significant pairwise differ- 
ences between sites across the three locations. The positive values 
of the neutrality tests, the presence of several common haplotypes, 
the smaller number of private ones, and the flat shape of the BSP, 
strongly indicate a long-term stability of the population size. This 
species may have persisted in two or more glacial refugia, followed 
by a secondary contact between them [5 7] . Although this ascidian 
has large dispersal capacity estimates [43] our results support 
larval retention, as has been previously suggested ([58], Figure 2), 
allowing genetic structuring in the range of 300 and 1500 km. 

Balanus perforatus, the species with the longest PLD estimates, 
lacked significant genetic structure across distances of up to 
1500 km. Neutrality tests indicated a very old expansion, starting 
around 60,000 ya. This pattern of persistence is similar to that 
shown by other widespread and abundant species inhabiting 
comparable habitats along the northeastern Pacific [59,60] . 

The genetic structure found among populations of Chiton 
olivaceus across the whole study area is congruent with a pattern 
of isolation by distance. This species showed the only estimate of 
expansion time clearly posterior to the LGM, indicating reduced 



persistence and recent colonization from glacial refugia [60] . This 
would also be in agreement with the low abundances and low 
molecular diversities detected with respect to the other species 
studied, and to other polyplacophora species [3]. 

As previously shown [61], our results support the weak 
correlation between PLD and genetic structure. The scale of 
larval dispersal is related to the complex interaction of larval life 
history traits, oceanographic regimes, habitat quality and distri- 
bution, and the variability of each of these factors through time. 
However, it is noteworthy that contrasting dispersal capacity 
patterns, as those of Hexaplex trunculus and Balanus perforatus, 
influence connectivity across the boundaries of biogeographic 
areas, as already detected in fishes [62]. The genetic structure 
observed is most likely due to the different effects of glacial and 
interglacial isolations and reconnections on their populations. 
Levels of gene flow are not the only factors determining present 
day genetic structure among populations; divergence times and 
population size can gready affect the observed genetic structure 
[63], however the lack of basic biological data for most of our focus 
species, does not allow to disentangle the role these different 
drivers play in the observed genetic structures. 

The two multivariate analyses allowed a general picture to be 
provided by merging in a common framework the results obtained 
for individual species. The statistical tests as well as the graphical 
representation provided a general picture of the phylogeographic 
structure in the Central Mediterranean, as well as the relative 
importance of the discontinuity areas detected, according to the 
number of species affected. This approach enables to summarize 
the available genetic information in a simple, straightforward 
result, which might be useful for decision-making or management 
purposes. Species with high haplotype and/or nucleotide diversi- 
ties could bias these analyses; however the multispecies approach 
counteracts this effect, so the results are not obscured by the 
structure of the most variable species. 

According to the results of the present study, the Tyrrhenian- 
Ionian biogeographical boundary acts as a phylogeographical 
barrier, as suggested by the genetic structure of two out of seven 
invertebrate species. The hierarchical sampling design adopted 
allowed the identification of an area of genetic discontinuity in the 
mid-Ionian Sea, reflected in the genetic structure of two species. 
The data retrieved from public databases reinforced the observed 
haplotype distributions, and gave further evidence on the 
variability of the effect of known barriers to gene flow on the 
genetic structure of different species. In the Mediterranean Sea, 
Pleistocene glaciations had a profound influence on the present 
day distribution of species and genetic diversity. The estimates of 
expansion and divergence allowed further understanding of the 
genetic structure found at present. Due to the ecological 
importance of the Central Mediterranean Sea as a hot spot of 
marine biodiversity [64], these results are of importance for the 
better understanding of the ecological functioning of its coastal 
ecosystem. General patterns pertinent to the management and 
conservation of this important ecosystem have been shown thanks 
to the comparative phylogeography approach used. Further 
research analyzing the genetic structure of other species, inclusive 
of additional Putative Boundaries, while also utilizing a larger set 
of molecular markers, would provide more insight on a multi- 
species, multi-scale picture of connectivity in the marine realm. 

Supporting Information 

Figure SI Mismatch distributions at each Location and 
for the whole data set for each of the seven species. Grey 
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lines represent the expected and black lines the observed 

distribution. 

(TIF) 

Figure S2 Bayesian Skyline Plots for the species show- 
ing signs of expansion according to mismatch distribu- 
tion or neutrality tests. Patella caerulea and Hexaplex trunculus 
plots are calculated for each haplogroup. When shown, the black 
square indicates the date of the Last Glacial Maximum (LGM). 
The most likely evolutionary model selected by jModelTest is 
indicated in each plot. 
(TIF) 

Table SI Specific primers. Names, sequence, amplified 
fragment length, and optimal annealing temperature for each of 
the seven pairs of primers used. 
(DOCX) 

Table S2 Previously published sequences. Accession 
numbers, geographical origin and identification with the haplo- 
types found in this study, as named in Figure 3, of sequences 
retrieved from Genbank. 
(DOCX) 
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